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ABSTRACT 

We investigate the effects of gravitational waves (GWs) from a simulated population of binary 
super-massive black holes (SMBHs) on pulsar timing array datasets. We construct a distribution 
describing the binary SMBH population from an existing semi-analytic galaxy formation model. Using 
realizations of the binary SMBH population generated from this distribution, we simulate pulsar timing 
datasets with GW-induced variations. We find that the statistics of these variations do not correspond 
to an isotropic, stochastic GW background. The "Hellings & Downs" correlations between simulated 
datasets for different pulsars are recovered on average, though the scatter of the correlation estimates 
is greater than expected for an isotropic, stochastic GW background. These results are attributable to 
the fact that just a few GW sources dominate the GW-induced variations in every Fourier frequency 
bin of a 5-year dataset. Current constraints on the amplitude of the GW signal from binary SMBHs 
will be biased. Individual binary systems are likely to be detectable in 5-year pulsar timing array 
datasets where the noise is dominated by GW-induced variations. Searches for GWs in pulsar timing 
array data therefore need to account for the effects of individual sources of GWs. 
Subject headings: black hole physics — galaxies: evolution — gravitational waves — methods: data 
analysis 



1. INTRODUCTION 

Detecting and performing science with gravitational 
waves (GWs) is currently a major goal of experi- 
mental astrophysi cs. Pulsar timin g arrays (PTAs, 
iHellings fc Downs 19831: iFoster fc Backer .199^ employ 
contemporaneous timing observations of millisecond ra- 
dio pulsars in order to search for the effects of GWs per- 
turbing the space-time metric along each pulsar-Earth 
line of sight. PTAs are complementary to other GW 
detection experiments, such as ground-based and space- 
based interferometers, in that they are sensitive to GWs 
in a different frequency band (currently ~ 5 — 100 nHz, 
i.e., periods of five years to a few months). In this fre- 
quency band, the most promising astrophysical sources 
of GWs are binary super-massive black holes (SMBHs). 

SMBHs, with masses in the range IO^Mq to \0^°Mq^ 
are of great importance to the evolution of galaxies. 
Feedback from active SMBHs is a key elemen t in shaping 
prope rties of the observed galaxy population (jGitti et al.l 
I2012f ). Models that trace the co-evolution of SMB Hs and 
their host galaxies fe.g.. [Pi Matteo et al.l [200 8^ postu- 
late that SMBHs grow primarily through accretion and 
coalescence with other SMBHs during active phases trig- 
gered by galaxy mergers. Presumably, there exists a large 
population of binary SMBHs at various stages of coales- 

^ Also at CSIRO Astronomy and Space Science, Australia 
Telescope National Facility, PO Box 76, Epping, NSW 1710, 
Australia 

^ E-mail address: v.vikram.ravi@gmail.com 
^ This mass range approximately corresponds to the cur- 
rent sample of S MBHs with dynamical mass measurements (cf. 
IMcConneli et al.ir2011 '). 



cence in the cores of galaxies that have recently merged 
with other galaxies. The final stages of SMBH-SMBH 
coalescence are driven by losses of energy and angular 
momentum to GWs, primarily emitted in the PTA fre- 
quency band. Various works have predicted the aver- 
age spectrum of the GW strain amplitude from the cos- 
mic population of binar y SMBHs ( Jaffc & Backer 200j; 
IWvithe fc L"oeb 2003; S esana et al.il2008D . Under the as- 
sumption that all binary systems are in circular orbits 
evolving o nly through G W emission, this spectrum takes 
the form (' Phinnevll200l 

/ic(/)=^lyr(^-^) ' , (1) 

where hdf) is the characteristic strain as a function of 
GW frequency, /, at the Earth, and Aiy^ is the charac- 
teristic strain at a frequency of /lyr = (lyear)"^. The 
characteristic strain is the GW strain per logarithmic 
frequency interval, and is given by 

hcif) = ^/JS^), (2) 

where Sh{f) is the one-sided strain power spectral den- 
sity (PSD). 

The combined GW signal from binary SMBHs is widely 
assumed to form an isotropic, stochastic GW back- 
ground. The value of Aiy^ is used to specify the am- 
plitude of the background. The m ost recent pred ictions 
for the value of Aiy^ were made bv lSesana et al.l p008i) . 
who found a likely range of 10~^^ to 2.5 x 10~^^. This 
range of predictions was derived by considering a variety 
of SMBH seeding, accretion and feedback scenarios, as 
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well as uncertainties in the galaxy merger rate and in the 
SMBH mass function. 

PTA projects are based on observations of periodic 
pulses of radio emission from pulsars using large radio 
telescopes. These observations lead to measurements of 
pulse times of arrival (ToAs) at the observatories. ToA 
measurements are typically conducted every few weeks 
over many years. Pulsar timing involves parameters of 
a physical model for the ToAs being fitted to the mea- 
sured To As, with the differences between the observed 
ToAs and the model predictions being the timing resid- 
uals. ToAs can be modeled using combinations of deter- 
ministic and stochastic processes. 

GWs incident on the Earth (and on the pulsars) cause 
shifts in the measured pulse frequencies of the pulsars 
([Sazhin 1978; Dctwcilcr 1979). For a pulsar, indexed by 
p, with intrinsic rotation frequency Vp, consider a GW- 
induced shift, Ai^(t, rp), to this frequency. This shift 
is a function both of time, t, and the Earth-pulsar di- 
rection vector, Fp. The resulting discrete time-series of 
GW-induced variations to the ToAs, 6f (the i subscript 
indicates that 5 f is sampled at times ti), is given by 
(|Detweilerill979ll 



-dt'. 



(3) 



For any GW signal, the expected values of zero-lag cross- 
correlations between the Sf time-series for different pul- 
sars can be specified. For any stochastic GW signal, 
the expected value of the normalized correlation for each 
pulsar pair is expressed in ter ms of the angular sepa- 
rations between t he pulsars as (jHellings fc Down^ll983l : 
IJenet et all 120051 ) 



3 , all, 
Ppq = -aloga- - + - + -Spg, 



(4) 



is the 



where a = 5(1 - cos 9pq), 9pq = cos ^ 
angular separation between pulsars p and q, and the 
Kroenecker delta, dpq, is unity ii p = q, and zero oth- 
erwise. This function is known as the Hellings & Downs 
curve. An isolated source of GWs will give rise to 
correlations that are different from the the Hellings & 
Downs curve. Measurements of correlations between pul- 
sar timing datasets that are attributable to the effects of 
GWs a rc necessary for the detection of GWs wit h PTAs 
(IJenet et al.l 120051 lYardlev et all 120111: [Demorest et al.l 
&)TW . 

The prospect of detecting or constraining the ampli- 
tude of a background of GWs from binary SMBHs has 
been the primary rationale for the development of the 
PTA concept. Various works have placed upper bounds 
on the value of Ai yr for a background with the charac- 
terist i c strain spectral form of Equation 1 (IJenet et al.l 
120061: Ivan Haasteren et al.l [20TTI : iDemorest et al.ll2012D 
The best published upper bound (jvan Haasteren et ahl 
1201 It ) finds that Aiy^ < 6 x lO'^^ with 95% confidence. 
All analysis methods developed to study the combined 
GW signal from binary SMBHs with PTAs assume that 
the 6f time-series for multiple pulsars can be described 
as a specific stochastic process. We describe the exact 
nature of this assumption in §2. 

In this paper, we elucidate the statistical nature of the 



ToA variations induced by GWs from binary SMBHs. 
We accomplish this by modeling the GW signal from the 
predicted population of binary SMBHs, and by simulat- 
ing realizations of Sf corresponding to realizations of the 
GW signal. This study is critical to the validity of inter- 
preting published upper limits on Ai yr as representative 
of limits on the mean characteristic strain spectrum of 
GWs predicted to arise from binary SMBHs. Our results 
are also important for the optimization of GW detection 
techniques with PTAs. In §3, we outline our method of 
simulating pulsar timing datasets including GWs from 
the predicted population of binary SMBHs. Our analy- 
sis and results are presented in §4 and §5, and we discuss 
the interpretation and implications of our results in §6. 
We present our conclusions in §7. 

Throughout this work, we assume a ACDM concor- 
dance cosmology based on a combi ned analy sis of the 
first-year WMAP data release (Soergel et al.l 12 0(3 3,') an d 
the 2dF Galaxy Redshift Survev (jColless et al.ll200l . 
with flM = 0.25, fib = 0.045, flA = 0.75, as = 0.9 
and Hq — 73kms~^ Mpc~^. Although these parameter 
values have since been superseded by more recent ob- 
servations, we adopt them in order to remain consistent 
with the model w e use for the binary SMBH population 
(|Guo et al.l 120111 ). A list of important symbols in this 
paper is shown in Table 1 , along with the sections of the 
text in which they are introduced. 

2. THE CURRENT MODEL FOR TOA VARIATIONS 
INDUCED BY GRAVITATIONAL WAVES FROM BINARY 
SMBHS 

ToA variations induced by GWs from binary SMBHs 
(^f) are commonly modeled among the PTA commu- 
nity as a wide-sense stationary random Gaussian process. 
This is based on the hypothesis that many GW sources 
forming a GW background contribute to the ToA vari- 
ations, resulting in a statistical process governed by the 
central limit theorem. While the nature of the random 
Gaussian mode l for 6f has been extensiv ely described 
elsewhere (e.g.. Ivan Haasteren et al.ll2009D . we summa- 
rize it here for completeness. 

The key property of a random Gaussian process is that 
a linear combination of samples will have a joint normal 
distribution function. Different samples need not be sta- 
tistically independent. The distribution of samples from 
a (zero-mean) random Gaussian process is characterized 
by the covariance matrix of the samples. Consider a vec- 
tor, Rp, containing n samples of 6f . That is. 
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Let Rq be another vector defined similarly to Rp, cor- 
responding to a pulsar q, containing n simultaneously- 
obtained samples of S'^. Under the random Gaussian as- 
sumption, the joint probability distribution of the sam- 
ples in Rp and Rq, which we denote Ppq, is given by 



p — 



1 



^(27r)"det(Cf 



(6) 



Here, Cpq is an n x rt matrix containing the covariances 
between the samples of Sf and (5|; that is, element ij of 
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Table 1 

List of symbols. 



Symbol Section 



Description 



5f 

Ppq 

SAf) 

ho 

? 
$ 

5g,fit(/) 
hc,Rtif) 

Wi 

SU) 

ppq 



1 
1 

2 
2 

3.1 

3.1 

3.1 

3.1 

3.2 

3.2 

4 

4 

4 

4 

5 



GW-induced ToA variation for pulsar p at time ti . 

Expected zero-lag normalized cross-correlation between <5^ and 5^ time-series. 

Expected PSD of 5^ time-series. 

Pcriodogram estimator of Sg{f) at frequency /fc. 

GW strain amplitude divided by frequency dependence. 

Binned distribution of binary SMBHs derived from realizations of the Millennium and Millcnnium-II coalescence lists. 
Average of 1000 realizations of <I>. 
Analytic fit to 

Sg{f) derived in terms of >l?fit. 

Expected GW characteristic strain spectrum derived in terms of $flf 

1 ns rms ToA variation at time ti . 

Sum of Wi and (5f . 

Expected PSD of time-series. 

Periodogram estimator of S{f) at frequency f^- 

Estimator of ppq. 



Cpq is given by the covariance between 6f and (5'. As the amphtudes set to be purely real with zero mean, variance 



Gaussian process is wide-sense stationary, each element 
ij of Cpq depends only on the time difference tij = \ti—tj\ 
between samples i and j for pulsar p and q respectively. 
Elements of Cpq are sampled from a covariance function, 
Cpq (r) , between the GW-induced ToA variations for pul- 
sars p and q. This covariance function is defined by the 
inverse Fourier transform of the one-sided PSD, Sg{f), 
of GW-induced ToA variations for a given pulsar: 



(7) 



Here, T denotes a complex Fourier transform, and r is a 
time-lag. The PSDs of the GW-induced To A variations 
for a ll pulsars are equivalent, and given by ()Jenet et alj 
l2006l) 

= 12^^. (B) 

for a GW signal with the expected characteristic strain 
spectrum hc{f)- 

The above discussion applies equivalently if pulsar p 
and pulsar q are the same pulsar, or if they are different 
pulsars. The Hellings & Downs factor ppq, defined in 
Equation 4, accounts for the correlation between GW- 
induced TOA variations for different pulsars. If hdf) has 
the form in Equation 1, we have Sg{f) oc f~^^^^. The 
GW-induced ToA variations for each pulsar will therefore 
be a "red" process. In this work, we only consider time- 
series with finite lengths Tobs- 

We are interested in comparing a new model for Sf with 
the random Gaussian model described in this section. To 
this end, we need to be able to simulate realizations of 
Sf as a random Gaussian process. Multiple PTA groups 
test their data analysis algorithms by sim ulating real- 
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izatio ns of Sf using the GWbkgrd plugin 
20091) to the tempo2 pulsar timing package 
2006) . While this plugin does not explicitly generate ran- 
dom G aussian r ealizati ons of Sf b y const ruction, we and 
others (jvan Haaster en et al.ll2011; Deni orest et al.|[2012D 
have checked that it approximates a random Gaussian 
process well. 

In the plugin GWbkgrd, a number of GW oscillators, 
Nt2, are simulated between GW frequencies /lo and /hi, 
with the normally distributed -|- and x GW polarization 



'ln(/hi//io) 



'T2 



T2 



(9) 



and frequency probability distribution, dP/df, given by 
dP_ l_ 

~W ~ HW/io) 



(10) 



ToA variations calculated for a given pulsar p at differ- 
ent times ti corresponding to GWs from each of these 
oscillators are summed to produce a realization of the Sf 
time-series. The frequency limits /lo and /hi are gener- 
ally chosen respectively to be much less than the and 
much greater than the Nyquist frequency corresponding 
to the minimum sampling interval. 

We make a distinction between the expected PSD of 
ToA variations induced by GWs from binary SMBHs, as 
defined in Equation 8, and estimates of this PSD based 
on realizations of the ToA variations. A commonly- 
used non-parametric, unbiased est imator of the P SD of 
a time-series is the periodogram (ISchusteriHSOl . The 

periodogram, Sg, of Sf is defined as 



CP 



|DFT[<5f]p, 



(11) 



where DFT denotes a discrete Fourier transform. We 
adopt the following standard definition for the DFT of n 
samples of Sf: 



DFT(/fc) 



n-l 
m— 



-i27rmk/n - 



(12) 



where i = \/— T in this case. The DFT is evaluated for 
frequencies 



/fc - (fc + 1) 



1 



To 



0< k < 



To 



bs 



bs 



2T 



(13) 



samp 



where Tsamp is the interval (assumed to be constant) be- 
tween samples of Sf. Throughout this work, we estimate 

the PSD, Sg{f), of realizations of Sf by evaluating S^. 
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3. SIMULATING PULSAR TOA VARIATIONS 
ACCOUNTING FOR BINARY SMBH POPULATION 
CHARACTERISTICS 

In this section, we describe a new method of simulating 
ToA variations caused by GWs from the predicted pop- 
ulation of binary SMBHs. Various works have presented 
models for the co smic demographics of binary SMBHs 
(jPotti et al.ll201^ . More recently, such efforts have been 
based on analytic prescriptions for baryon physics ap- 
plied to dark matter halo m erger tree catalogues from 
N-body simulations ()Guo et al. 2011,, hereafter Gil). In 
particular, the Gil prescriptions w ere applied to merger 
trees from both the Millennium (iSpringel et al.l 20051 ) 
and the Millennium-II (|Bovlan-Kolchin et al.ll2009 ) sim- 
ulations. The Millennium and Millennium-II simulations 
follow the evolution of dark matter structures, using the 
same physical prescriptions and number of particles. The 
Millennium-II simulation was however carried out in a 
comoving cubic volume with one-fifth the side length 
as that of the Millennium simulation, with the aim of 
resolving smaller-scale dark matter structures than the 
Millennium simulation^ Together, these simulations re- 
solve dark matter halos corresponding to the observed 
galaxy population, from dwarf galaxies to the largest- 
mass early-type galaxies. 

Th e Gil model is the latest in a series ( Sprinac l et al.l 
[200l ICroton et al.1 [20061 IDe Lucia fc Blaizota i2007l ) of 
semi-analytic prescriptions applied to the Millennium 
simulations. A host of observables of galaxies at low red- 
shifts are reproduced, along with the redshift-evolution 
of the quasar population and star formation. Of most 
relevance here is that the model also traces the SMBH 
population, reproducing the 2 = SMBH-galaxy scaling 
relations in their slopes, normalization and scat ters, as 
well a s the inferred SMBH mass function (Marull Tet al.l 
[2008h . 

We base our description of the binary SMBH popu- 
lation emitting GWs on the prediction for the SMBH- 
SMBH coalescence rate from the Gil model. We fit an 
analytic function to the distribution of binary SMBHs, 
and randomly draw GW sources from this distribution to 
produce realizations of the GW sky corresponding to bi- 
nary SMBHs. We then add the effect of each GW source 
to simulated pulsar ToA datasets in order to analyze the 
GW-induced ToA variations. 

This work is different from previous attempts to model 
the GW signal from bina ry SMBHs. Initial attempts 
(e.g.. lJaffe fc Backeil200l to predict the mean GW char- 
acteristic strain spectrum from binary SMBHs used em- 
pirical determinations of the galaxy me r ger ra te and the 
SMBH mass function. iWvithe fc Loebl (|2003D predicted 
the GW spectrum by analytically following the dark mat- 
ter halo merger hierarchy in an extended Press-Schechter 
framework, and by deriving the SMBH coalescence rate 
by re lating the SMBH ma sses to the halo circular veloc- 
ities. iSesana et al.l ()2008[ ) considered the possible range 
of predictions of the characteristic strain spectrum, us- 
ing Monte Carlo realizations of dark matter halo merger 
trees and various prescriptions for SMBH growth. 

The key difference between the present work and pre- 
vious calculations of the GW signal from binary SMBHs 

* The Millennium-II simulation however does not reproduce 
larger-scale structures as well as the Millennium simulation. 



is that we are chiefly concerned with the statistics of 
5^. Our approach to modeling the bi nary SMBH popu- 
lation is similar to ISesana et al.l (|2009[ ) in our use of mock 
galaxy catalogues derived from analytic prescriptions ap- 
plied to the Millennium simulations. However, whereas 
Sesana et al.l ([2009) modeled the SMBH population by 
using empirical SMBH-galaxy scaling relations combined 
with (earlier) mock catalogues, we utilize SMBHs mod- 
eled by Gil in a self-consistent framework which repro- 
duces the relevant observables. 



3.1. Modeling the distribution of binary SMBHs 

As in the previous works discussed above, we consider 
all binary SMBHs to be in circular orbits, and use expres- 
sions f or the resulting GW emission presented bv lThorni 
(flOSl . We briefly discuss the assumption of circular or- 
bits in §6.3. The strain amplitude, hs{f), at frequency / 
of GWs from a circular binary, averaged over all orienta- 
tions and polarizations, is given by: 



hsif) 



d^D{z) 



(7r/(l + z))2/3, (14) 



where G is the universal gravitational constant, Mc = 
(MiM2)^/^(Mi-hM2)"^/^ is the chirp mass of the binary, 
c is the vacuum speed of light, z is the redshift, and D{z) 
is the comoving distance to the binary. The evolution 
of the received GW frequency with observed time, t, is 
determined by 



dt 



V«/3/"/'(GMc(l + z))'^/^ (15) 



The rest-frame binary orbital f requency is giv en by /f, = 
i/(l + z). We assume, after IHughesI (|2002D , that the 
maximum received frequency, /max, is attained at a bi- 
nary separation of three Schwarzschild radii, correspond- 
ing to the last stable orbit: 



/n 



12\/37r(l + z)GM, 



(16) 



c 



assuming a mass ratio of unity. 

The mock galaxy catalogues resulting 
from the Gil model are available onlinCi 
(|Lemson fc Virgo ConsortiumI 120061) . The halo merger 
trees from the Millennium and Millennium-II simulations 
were evaluated at 60 logarithmically-spaced redshift 
"snapshots" between z = Q and 2; = 20. We obtained 
the lists of SMBH-SMBH coalescence events within the 
comoving volume of each simulation by querying the 
online database. Rcdshifts at the (non-logarithmic) 
midpoints between the redshift snapshots were assigned 
to each event. We used these lists to fill bins of a 
distribution, $, of the number, N , of observable binary 
SMBHs per unit comoving volume per solid angle on 
the sky, given by 



dN d^Vc dz dt 

$ = An 

dhci dfldz dt df 

^ http:/ /www. mpa-garching.mpg.de/millennium/ 



(17) 
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where 

(18) 

and 47r^j^ is the sky-integrated comoving volume 
shell between redshifts z and dz. Also, % = 

' at 

iJo\/f^Af(l + ^Y" + ^A, and the derivative ^ was ob- 
tained from Equation 15. $ is the predicted distribu- 
tion of binary SMBHs in Hq (which corresponds to the 
frequency-independent GW "power" , or squared strain 
amplitude) and the observed GW frequency. 

For chirp masses below 10 ''M0, the limited capability 
of the Millennium simulation to resolve low-mass halos 
caused an under-prediction of the chirp mass function 
as compared to the Millennium-II simulation. In order 
to ensure a complete chirp mass function, we included 
binary SMBHs with Mc > lO^M© from the Millennium 
list of coalescence events, and binaries with lO^Af© < 
Mc < 10"^ Mq from the Millennium-II list. 

Some degrees of randomization in the coalescence lists 
were possible. First, in cases where more than two 
SMBHs coalesced to form a single SMBH between red- 
shift snapshots, the merger order was not specified. In 
these instances we randomized over merger order. Sec- 
ond, a spherical comoving volume shell between any pair 
of redshifts less than ^ 0.09 could be contained within 
the simulation volume. Some Millennium redshift snap- 
shots exist at z < 0.09, and the comoving volume shells 
between redshifts corresponding to these snapshots en- 
close some SMBH-SMBH coalescence events in the Gil 
model. An observer located at the center of the Mil- 
lennium simulation volume would observe only a frac- 
tion of the total list of events in the Gil model at 
z < 0.09, and an observer located elsewhere in the vol- 
ume would observe a different selection of events. This 
is not the case for the Millennium-II simulation, where 
the volume was too small to enclose any comoving vol- 
ume shells between redshift snapshots. For each real- 
ization of the Millennium (but not the Millennium-II) 
coalescence list, we therefore specified randomly-placed 
spherical shells within the simulation box to select bi- 
nary SMBHs at these redshifts. For coalescence events 
at 0.09 < z < 0.19, the corresponding comoving volume 
shells between redshift snapshots were smaller than the 
Millennium simulation volume, though not enclosed by 
it. For these coalescence events, we randomly included 
binaries in the Millennium coalescence list according to 
probabilities given by the ratios between the volumes of 
the comoving shells and the Millennium simulation vol- 
ume. We used 1000 realizations of the Millennium and 
Millennium-II coalescence lists to form realizations of the 
binary SMBH distribution $. 

In generating realizations of the distribution $, we 
assumed that every SMBH-SMBH coalescence event in 
the Gil model catalogues was the result of a binary 
SMBH system that had decayed through GW emission. 
The Gil model included the assumption that, upon the 
merger of two galaxies with central SMBHs, the SMBHs 
coalesced in every case, before accretion onto the newly- 
formed SMBH0 We note that SMBHs with masses as 

^ There are various mechanisms by which extreme mass ratio 



low as IO^Mq were present in the Gil model catalogues, 
but were not included in the $ distributions. We veri- 
fied that relaxing the lower cutoff on the SMBH masses 
in the $ distributions from IO^Mq to lO^M© did not 
significantly modify the total signal. 

The 1000 realizations of $ were averaged to form a 
distribution We fitted $ with an analytic function 
which could be used to generate random realizations of 
the observable binary SMBH population. We did not 
use realizations of $ as realizations of the binary SMBH 
population because the ^-distributions were binned for 
computational purposes. A four-parameter function, 

with free parameters n, ph, a and /?, was found to fit $ 
well. We performed the fit on the logarithm of the data 
to approximate linearity in the fitting procedure. The 
best-fit parameters are given in Table 2. The frequency- 
exponent was held fixed at —11/3, as predicted by Equa- 
tions 15 and 17. 

3.2. Realizations of pulsar ToAs with GW-induced 
variations 

For a set of binary SMBHs drawn from the distribution 
<l>fit , we simulated a corresponding time-series by sum- 
ming the contributions from each individual binary. De- 
tails of the meth od used to ca l culate these contributions 
are presented in ' Hobbs et al.l (|2009f ) . For each binary, 
we randomized over the right ascension and declination, 
the orbital inclination angle, the orientation of the line 
of nodes, and the orbital phase angle at the line of nodes. 
A new pubhcly-available tempo2 plugin, "addAllSMB- 
HBs" , was written to perform this simulation. 

For most of the present work, we did not use tempo2 
to fit timing model parameters. Instead, we made use of 
the tools available for spectral analysis of timing resid- 
uals. In our simulations, the "timing residuals" corre- 
sponded exactly to 6f given the absence of timing model 
fitting. 

We consider it important to emphasize the distinction 
between the 5f time-series and the timing residuals re- 
sulting from analyses of observed ToA datasets. Consider 
a set of observed ToAs that exactly match a particular 
timing model, except for the addition of GW-induced 
variations (a time-series) . Given that an observer does 
not actually possess any prior knowledge of the timing 
model parameters, the observer will fit the model pa- 
rameters to the ToAs. The resulting timing residuals 
will not be equivalent to . This is because the vari- 
ations in the ToAs can alter the apparent pulsar timing 
parameters. For example, the presence of a Sf time-series 
consisting of a sinusoidal signal with a period of one year 
will alter the apparent pulsar position. 

In order to investigate the statistics of given our 
model for the binary SMBH population, we first used 
TEMP02 to generate 500 ToAs spanning 5 years exactly 
corresponding to the P SR J0437— 4715 timing model 
([Manchester et aH l2012| ) . We then added realizations 
of the 6f time-series evaluated at the observed ToAs to 

binary SMBH system s and triple or hig her-order systems can avoid 
coalescence (e.g.. .Volonteri et al.|[2003 ). 
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Table 2 

Best- fit parameter values of ^at- 



Parameter 


Value 


n 


2087 ± 365 


Ph 


4.878x10-23 ±4.45 x IQ-^^ 


a 


-1.72249±0.00064 


/3 


-0.3473±0.0046 



these datasets, that is, with Tobs = 5 yr and Tsamp = 
0.01 yr. The pulsar distance was set to Ikpc, and the 
position was held fixed for all simulations. As the binary 
SMBHs used to produce realizations of Sf had random- 
ized positions and orientations, allowing the pulsar posi- 
tion to vary between realizations would not alter our re- 
sults. The results presented in this paper are not depen- 
dent on the timing model used or on the pulsar distance 
from the Earth. We also added Gaussian white noise 
variations with 1 ns rms to the ToAs. This white noise 
component is miich smaller than is usually observed in 
ToA datasets, but was necessary to smooth over machine 
precision errors. 

We included GW sources between 10~^ Hz and 10~^ Hz 
in our simulations of Sf. The lower frequency cutoff was 
chosen to be less than one fifth of /o = (5 years) The 
upper frequency cutoff was chosen to be greater than 
,/249 = (0.02 years) We assumed, after previous works, 
that all GW sources between these frequency cutoffs are 
non-evolving over a 5-year timespan, i.e., they do not 
evolve in frequency by more than (5years)~^. The up- 
per bound on the /iQ-values of sources, ho, max(/), was set 
by the last stable orbit of binary SMBHs. We identified 
this bound by fitting a power law to the high-/io edge 
of the $ distribution. The lower bound on ho, /iq, min, 
was set by the lowest non-zero /lo-value in This value 
corresponds to a binary SMBH containing two lO^M© 
components at z « 6. The distribution included more 
than 6.5 x 10^® GW sources within this ho — f domain; 
the vast computational cost involved makes it impossi- 
ble to simulate Sf with this many sources. Fortunately, 
the shape of the $fit distribution was such that, at a 
given frequency, the highest- ft,o sources contributed most 
to Sg{f) (we return to this point below), defined in terms 
of $f;t (using Equation 8) as: 



^g,fit(/) = 



'lO,max(/) 



127r2/' 



/■/iO, me 

ho^ min 



^^thl{f)dho, (20) 



The average characteristic strain spectrum derived from 
$fit is 



hoMf) = f 




ho, max(/) 



1/2 



$fit/i'(/)d/io 



(21) 



We found a function, ho{f ), such that 
.(/) 

0.95g,fit(/) 



rn-o, I 



dho = 5'g,fit(/)- 



hoif) 12^V^ 

(22) 

Thus, the GW sources in the domain /io(/) < ho < 
^o,max(/) contribute, on average, 90% of Sg{f) at every 
frequency. Between 10~® Hz and 10~^ Hz, this amounted 
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Figure 1. Illustration of the ho — f domain constraints on the 
distribution '^tit- The upper and lower dashed lines represent 
/io,max(/) and Ao.min, as labelled, and the solid curve represents 
ho{f). The shaded region is the "90% domain" from which binary 
SMBHs contributing, on average, 90% of the ToA variation PSD 
at every frequency were drawn. 

to ~ 4.5 X 10^ sources. We refer to this ho ~ f do- 
main as the "90% domain". The 90% domain, along 

with /io,max(/), hoif) and fto.min, is shown in Figure 1. 

We approximated the total number of sources (6.5 x 
10^^) in the ho — f domain between ho. min and /lo.max 
and between 10~^ Hz and 10~® Hz as constant. For a 
given realization of the actual number of sources in 
the 90% domain is governed by binomial statistics. We 
therefore drew a (binomial-) random number of sources 
from the 90% domain, and added contributions from 
each of them to each realization of (5f . We assumed 
that the sources remaining in the distribution with 

^0, min < ho < ho{f), Contributing on average 10% to 
Sg{f) at every frequency, resulted in a stochastic contri- 
bution to 6^ governed by the central limit theorem. We 
therefore simulated them as described in §2 using the 
method of simulating ToA variations corresponding to 
a GW background implemented in the tempo2 plugin 
GWbkgrd. We simulated Nt2 = 5 x 10'* sources be- 
tween 10^^ Hz and 10~^ Hz using the tempo2 method, 
with the characteristic strain spectrum given by hc{f) = 

^127r^/^(S'g,fit(/) — Sg,fit{f)). For each realization of 

Sf, we added contributions from the ^ 4.5 x 10® GW 
sources drawn from the 90% domain of the ^fjt distribu- 
tion, and from the 5 x 10'* GW sources corresponding to 
the remaining (on average) 10% of Sg{f) drawn using the 
TEMP02 method. Shifts in the measured pulse frequen- 
cies caused by metric perturbations at both the Earth 
and the pulsar (i.e., the "Earth term" and the "pulsar 
term") were included in our simulations. 

The top panel of Figure 2 shows /ic,fit(./) in the 
< fc < 100 spectral bins. We also show a charac- 
teristic strain spectrum of the form in Equation 1 with 
^lyr = 7.8 X IQ-^^. This value of Aiy^ can be taken to 
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Figure 2. Top: The solid curve shows the mean characteristic 
strain spectrum, /ic.fiti derived from the distribution <l>fit in Equa- 
tion 21. The dashed line shows a representative spectrum of the 



form in Equation 1, with Ai^ 



7.8 X 10" 



The values of both 



traces are equivalent at the lowest frequency. The dotted line shows 
a spectrum of the form in Equation 1 with Aiyr = 6 X 10~^^, cor- 
responding to t he most recently publishe d 95% confidence upper 
bound on ^lyr Ivan Haasteren et al.l2011l) . Bottom: The numbers 
of GW sources that contribute 50% (dashed line) and 90% (solid 
line) of 5g_flt(/)- The numbers are integrated over frequency bins 
of width (5 years)^-*^ Hz. 



be the prediction from the Gil model, as it corresponds 
to the characteristic strain in the lowest (k = 0) spectral 
bin of a 5-year dataset. This particular curvat i ire in the 
/ic, fit curve, also predicted bv lWvithe fc Loebl (|2003[ ). is 
caused by the frequency-dependence of /io,max(/), which 
represents the bound beyond which binary SMBHs have 
crossed the last stable orbit. This curvature is caused 
by different p roces ses from the curvature reported by 
iSesana et all ()2008[ ). We show the mean characteris- 
tic strain spectrum for GWs fro m all binary SMBHs 
in the predicted distribution $fit. ISesana et al.l ()2008f ) 
fitted a broken power-law to realizations of the char- 
acteristic strain spectrum, accounting for various ran- 
domizations o ver th e source population. In particular, 
ISesana et al.l ()2008D randomized over the existence of 
"fractional" sources in every frequency bin of a fiducial 
dataset, and also excluded the strongest single source in 
every frequency bin in an attempt to isolate the back- 



ground signal. The smaller number of sources per unit 
frequency at higher GW frequencies, combined with the 
greater contributions to the signal from the strongest sin- 
gle sources in frequency bins at higher frequencies, both 
resulted in the curved c harac teristic strain spectra pre- 
sented by S esana et al.l (j2008[) . 

The bottom panel of Figure 2 shows the mean num- 
bers of highest-/io sources that contribute 90% and 50% 
of "Sg, fit(/) in these frequency bins. A small num- 
ber of sources contribute a large fraction of Sg^mif) 
at every frequency. In the fc = frequency bin, the 
^ 3 X 10"* highest-ft-o sources contribute on average 90% 
of Sg^titif), and only 30 sources on average contribute 
50% of S'g,fit(/)- At frequencies / > 1.5 x 10"'^ Hz, the 
strongest source, on average, contributes more than 90% 
of Sg^fitif) in each frequency bin. This is a consequence 
of the shallow power-law nature of the ft,o-distribution of 
the GW sources in the 3>fit distribution. 

In this work, we compare the Millennium-based sim- 
ulations of Sf with simulations of 6f created using the 
TEMP02 method described in §2. To this end, we sim- 
ulated ToAs as before, but added realizations of 6f cor- 
responding to 5 X IC* oscillators simulated using the 
TEMP02 plugin GWbkgrd. These oscillators were sim- 
ulated as described in §2 with a mean characteristic 
strain spectrum given by /ic, fit(/)- We refer to this latter 
method of simulating Sf as Case H09, after iHobbs et al.l 
(pool . Simulations of using GW sources drawn from 
$fit will be referred to as Case R12 after the present 
work. 

4. FOURIER-SPECTRAL ANALYSIS, AND RESULTS 
In this section, we consider the differences between the 
cases in the distributions of the periodograms, S*^, evalu- 
ated for realizations of Sf for a single pulsar. This is mo- 
tivated by the results in Figure 2, in particular that the 
number of GW sources per spectral bin that contribute 
90% of S'g, fit(/) varies significantly with frequency. The 
Case R12 simulations are intended to more accurately 
represent the effects of GWs from binary SMBHs on ToA 
datasets than the Case H09 simulations. Example real- 
izations of 5-year 6f time-series in both cases are shown 
in Figure 3. The time-series appear quite similar: real- 
izations in both cases are dominated by low-frequency 
components. Values of up to 1 /is are also present in one 
realization. 

Instead of directly measuring S^, the added white noise 
component in the simulated ToAs required us to analyze 
the periodograms of a time-series, Df , given by 



(23) 



where Wi is a time-series of Gaussian white 1 ns rms 
ToA variations as discussed above. The PSD of , 
Sg,fit{f), is significantly red, with a spectral index of 
— 13/3 (see Equations 1, 20), and is expected to domi- 
nate the PSD of Wi at low frequencies. We u sed the gen- 
eralized least-squares algorithm described in iColes et al.l 

()2011| ) to measure the periodograms, ipk of realizations 
of Df. This method requires knowledge of the auto- 
covariance function of Df , which we obtained as in Equa- 
tion 7 using the inverse DFT of the known PSD of Df , 
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Days 

Figure 3. Example realizations of 5? in Case H09 (thick grey 
lines) and Case R12 (thin black lines). 



S{f), given by 

sif) = SsMf) 



2(lns)2 



250/(5 years) 



(24) 



In the following, we consider the distributions of ip^ in 
the lower spectral bins, where S{f) w S'g,fit(/), to be 

approximately equivalent to the distributions of S^.. 
We produced 1000 realizations of Df in Case R12 and 

in Case H09, and measured f/'^ for each realization. In 

the top panel of Figure 4, we show the averages of 
measured from each of the Case R12 and Case II09 re- 
alizations, along with the expected PSDs of 6f and Wj. 

Arbitrarily chosen single realizations of t/j^ in each case 
are also shown. The means of the pcriodograms in both 
cases arc clearly the same, and equivalent to the pre- 
dicted PSD, S{f), given in Equation 24. Though this is 
as expected, it is both a check of the simulations of D^, 
and a demonstration of the ability of our PSD estimation 
method to measure steep red spectra without bias. 

In contrast, the distributions of tp^ in the frequency 
bins where S{f) « 5g,fit(/) are different between the 

cases. The single realizations of ijj^ in each case shown 
in the top panel of Figure 4 begin to hint at these differ- 
ences. In most spectral bins, the Case R12 periodogram 
is below the Case H09 periodogram. That this is a gen- 
uine trend is confirmed in the bottom panel of Figure 4. 
Here, we depict various "percentile" pcriodograms of the 

distributions of in each of Case R12 and Case HOQ. 
The percentile pcriodograms may be interpreted as con- 
tours of equivalent percentiles of the periodogram dis- 
tributions in different spectral bins. For example, the 
"50%" percentile periodogram links the 50th percentile 
points of the distributions of periodogram values in each 
spectral bin. Below the 95th percentile, all Case R12 per- 
centile periodograms lie below Case H09 percentile peri- 
odograms. This implies that in most spectral bins, most 
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Figure 4. Top: The mean estimates (ipl.) of the PSD of the 
simulated ToA variations (Df ) in Case R12 (thin solid bla<;k line) 
and in Case HOQ (thick solid grey line). The predicted PSDs of 

('S'g. fit (/)) ^nd arc shown as sloped and horizontal dashed 
lines respectively. Randondy-choseu single measurements of in 
Case R12 and Case HOQ are also shown, scaled down by a factor 
of 10, as black and grey dotted lines respectively. Bottom: The 
thin black and thick grey curves depict "percentile periodograms" 
of the distributions of Case R12 and Case HOQ measurements of 

respectively. The 5th, 25th, 50th, Q5th percentiles are shown 
as labelled, along with the maximum values of the periodograms 
in each spectral bin (labelled "max"). The vertical dashed lines 
indicate the k = and fc = 10 spectral bins, with frequencies given 
by (fe + 1)(5 years)-! Hz. 

measurements of a periodogram in Case R12 will, like the 
individual ones shown in the top panel of Figure 4, be 
below most Case H09 periodograms. However, the 95th 
percentile periodograms are essentially equivalent, and 
the maximum value Case R12 periodogram is well above 
the maximum value Case H09 periodogram. These im- 
plied "tails" at high values in the Case R12 periodogram 
distributions in each spectral bin are highlighted in Fig- 
ure 5, which depicts the distributions of the tp^ in the 
spectral bins indicated by the vertical lines in the bot- 
tom panel of Figure 4. The distributions are shown as 
the fractions of Case R12 and Case H09 periodograms at 
or above a given value. 

Figure 5 also shows that the Case R12 periodogram 
distribution in the fc = 10 spectral bin has a longer tail 
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Figure 5. The distributions 1000 measurements of in Case 
R12 (thin black lines) and in Case H09 (thick grey Hnes) in the 
A; = (top) and fc = 10 (bottom) spectral bins. The distributions 
are shown as the fractions of realizations at or above a given value. 
The domains of both plots indicate the maxima and minima of the 
distributions. 

relative to the Case H09 distribution, as compared to the 
fc = spectral bin. This effect is also evident in the bot- 
tom panel of Figure 4, in that the fractional differences 
between the percentile periodograms are greater at the 
upper end of the GW-dominated frequency regime. This 
is consistent with the number of GW sources per spec- 
tral bin included in the Case R12 simulations going down 
with increasing frequency, as shown in Figure 2. 

In summary, approximating ■0^ with as discussed 
above, we find that: 

• In most spectral bins, most realizations of in 

Case R12 will be below most realizations of in 
Case H09. 

• The maximum possible values of in Case R12 
will be higher than the maximum possible values 
of Si in Case H09. 

5. ESTIMATES OF CORRELATIONS BETWEEN 
GW-INDUCED TOA VARIATION TIME-SERIES 



iHellings fc Dowiisl (|1983[ ) showed that the average val- 
ues of correlations between for different pulsars, for 
a stochastic GW signal, will always be given by the 
Hellings fc Downs curve (Equation 4). The distribu- 
tions of individual estimates of these correlations, much 
like the distributions of the PSD estimator considered 
above, will however depend on the nature of the GW sig- 
nal. Here, we characterize the distributions of estimates 
of these correlations for multiple pulsar pairs in each case 
discussed above. We simulated 100 realizations of each 
of Case II09 and Case R12 ToAs as described in §3.2 for 
pulsars at the positions of each of the 20 pulsars timed 
by the Parke s PTA, using th e timing models specific to 
each pulsar (jManchester et al. 2012). For each realiza- 
tion, the same set of GW sources was used to simulate 
GW-induced ToA variations for each pulsar. The pulsar 
distances were set arbitrarily between Ikpc and 20kpc. 

We estimated the correlations between time-series i5f 
and , ppq, for each pulsar pair pq in each realization 
of Case R12 and Case H09 To As. No autocorrelations 
were estimated. A frequenc y-domain estimation tech- 
nique, based on the method of lYardlev et al.l (|2011| ). was 
used. This technique will be detailed elsewhere (Hobbs 
et al., in preparation). We refer to our estimates of ppq 
as Ppq. These estimates were performed using time- 
series, rather than 5f time-series (see Equation 23), and 
the estimation technique was optimized using the ex- 
pected PSD of given in Equation 24. The technique 
removes best-fit linear and quadratic terms from each 
time-series using the standard tempo2 least-squares 
fitting algorithm. This mimics the effect of fitting pulse 
frequency and frequency-derivative terms to the simu- 
lated ToAs and then analyzing the timing residuals. 

Following the removal of linear and quadratic terms 
from one of the 100 Case R12 realizations of D^, a single 
GW source was found to dominate the residual time- 
series. We show the corresponding Df time-series for two 
of the 20 simulated pulsars for this realization in Figure 6. 
The left panel of this Figure shows the large sinusoidal 
oscillations induced by the source, and the right panel 
shows example Case R12 realizations of that are not 
dominated by an individual source. It is possible that 
an individual GW source with a period greater than the 
5-year dataspan could dominate the realizations of Df 
in the right-hand panel of Figure 6. The ToA variations 
induced by such a source would however be absorbed in 
the removal of the linear and quadratic terms from the 
Df time-series. 

We averaged all measurements of ppq for each pulsar 
pair pq from the Case R12 realizations, besides the one 
clearly dominated by an individual source. The realiza- 
tion dominated by an individual source added a large 
amount of scatter to the average Case R12 correlations, 
and was left out of the average to enable a better compar- 
ison between the cases. We also averaged the Case H09 
measurements of ppq for each pair pq in 99 arbitrarily- 
chosen realizations. The average measurements of ppq are 
shown for both cases in Figure 7. The functional form 
of the Hellings fc Downs curve is recovered in both Case 
R12 and Case H09. However, the Case R12 estimates are 
significantly more scattered about the expected values of 
the correlations than the Case H09 estimates. 

The increased scatter in the Case R12 correlations with 
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Figure 6. Two examples of simulated realizations of £>? for two pulsars: PSR J1600— 3053 and PSR J1909— 3744 (see text for details). 
The realization in the left panel is affected by a strong individual GW source, whereas the realizations in the right panel is not. The lower 
plots show the time-series from the corresponding upper plots with linear and quadratic terms removed. 
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Figure 7. The estimated correlations between GW-induced ToA variations for simulated pulsars at the positions of the 20 Paikes PTA 
pulsars in Case H09 (left) and in Case R12 (right), plotted against the angular separations on the sky between each pair of pulsars. Each 
point represents an average over 99 realizations; in Case R12, one realization including an extremely strong individual source was not 
included in the average. Linear and quadratic terms were removed from each ToA variation time-series. The solid curve is the expected 
Hellings & Downs curve given in Equation 4. As no autocorrelations were present, the maximum value of the Hellings & Downs curve is 
0.5. 



respect to the Case H09 correlations in Figure 7 is caused 
by outlying estimates in only a few realizations of ToAs. 
This is shown in Figure 8, where we display the his- 
tograms of the ppq measurements between simulated ToA 
datasets for PSR J0437-4715 and PSR J0613-0200 in 
each case. Correlation estimates \ppq\ > 1 were possi- 
ble because we normalized the estimated covariances be- 
tween time-series using the expected cross-PSD be- 
tween the time-series. While most measurements in both 
cases are concentrated about the expected value of ppq, a 
few Case R12 measurements are significantly displaced. 
This is consistent with the results of §4. We also stress 
that the large scatter of the estimator common to both 
cases is expected, and intrinsic to the GW signal. 



6. DISCUSSION 

We have shown that the ToA variations induced by 
GWs from the predicted binary SMBH population are 
not consistent with the model described in §2. The ran- 
dom Gaussian model for 5^ described in §2 and approxi- 
mated in Case H09 is reasonable given that a large num- 
ber of GW sources are expected to contribute to the 
GW-induced ToA variations. That is, the values of 6^ 
at all times U are the sums of many random variables. 
An argument based on the classical central limit theorem 
would suggest that 5f would then be Gaussian random 
at every time ti. It is apparent, however, that such a 
central limit theorem-based argument does not apply to 
the Case R12 realizations of 6f. This is because of the 
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Figure 8. Histograms showing the distributions of measurements 
of the correlation estimator ppq for 100 simulated ToA datasets for 
PSR J0437-4715 and PSR J0613-0200, in Case R12 (top) and 
Case H09 (bottom). See the text for details of the simulations. 
The Case R12 realization that included an extremely strong single 
GW source, as discussed in the text, resulted in a measurement 
of Ppq = —39.74; this measurement is not shown in the Case R12 
histogram. The vertical dashed line in each panel indicates the 
mean values of all 100 estimated correlations in each case, and the 
vertical dotted line indicates the expected value of the correlation, 
Ppq, for an angular separation of 6pq = 49.8°. 

nature of the GW sources contributing to 5f in Case R12. 

In Case R12, a few sources contribute most of the PSD 
of Sf at every frequency, as shown in Figure 2. These 
sources are rare, because they are found at the high-Zip 
tail of the $fit source distribution. The estimators we 
consider in this work, and ppq, are dominated in Case 
R12 by a few GW sources that need not occur in every 
realization of the 5^ time-series. This is why the distribu- 
tions of these estimators are different between the cases. 
The quantities that we estimate, Sg{f) and ppq, are used 
to define the covariance matrix of the GW-induced ToA 
variations (see Equation 7). We have therefore shown 
that the ToA variations induced by GWs from binary 
SMBHs are dominated by the effects of a few strong, 
rare sources and cannot be accurately modeled using the 
random Gaussian process discussed in §2. 

6.1. Implications of our results for experiments focused 
on a GW background 

Current PTA data analysis techniques use assump- 
tions about the statistics of to attempt to estimate 



or constrain the amplitude of the characteristic strain 
spectrum of GWs from binary SMBHs. We consider the 
implications of our results for a selection of techniques 
in turn. We assume, in this discussion, that our results 
for the statistics of GW-induced ToA variations would 
apply even if the normalization of the GW characteristic 
strain spectrum /ic,fit(/), which we refer to as the GW 
amplitude^ were scaled up or down. Such a scaling could 
occur, for example, under different scenarios for whether 
coalescing SMBHs accrete gas before or after coalescence 
(|Sesana et a l, 2008). 
We summarize a few key techniques here: 

• iJenet et al.l ([20051) describe a statistic which mea- 
sures the degree of correlation between estimates 
of ppq from ToA data, and the expected functional 
form of Ppq. The expected detection significance, 
which is estimated under the assumption that the 
GW-induced ToA variations are Gaussian random, 
saturates at high values of the GW amplitude once 
the variance of the statistic is dominated by the 
stochasticity of the GW signal (see our Figure 8 
and related discussion). 

• IJenet et al.l (|2006D constrain the amplitude of the 
GW characteristic strain spectrum from binary 
SMBHs by estimating the maximum possible GW 
signal present in measured data, under the assump- 
tion that the data could be modeled using a white 
noise process and GW-induced ToA variations. A 
statistic that estimates the GW background am- 
plitude from individual pulsars was measured, and 
compared to the simulated distributions of the 
statistic for different GW amplitudes. The simu- 
lated statistic distributions were created from sim- 
ulated ToA datasets with GW-induced ToA varia- 
tions included using the tempo2 plugin GWbkgrd. 

• Ivan Haasteren et al.l ()2009D present a Bayesian 
parameter estimation method for the GW 

characteristic stra in spectral amplitude^ 

Ivan Haasteren et al.l ()2011[ ) used this method 
to constrain the GW amplitude. This method 
requires an evaluation of the likelihood of the 
parameters used to model the ToA datasets, which 
include the GW amplitude. The likelihood is 
the probability distribution of the data given the 
model parameters. The GW amplitude is used 
to calculate the covariance matrix, Cpq, between 
the GW-induced ToA variations for pulsars p and 
q (see Equations 7, 8). This covariance matrix 
in turn is used to define the PTA likelihood, 
assuming that the GW-induced ToA variations 
can be modeled as a random Gaussian process. 

• iDemorest et al.l (1201 21) use a PTA likelihood sim- 
ilar to the work of Ivan Haasteren et aD (|2009[ ) to 
constrain the GW amplitude by evaluating the dis- 
tribution of a maximum likelihood estimator for 

^ We make a distinction between this amplitude and the Aiyr 
parameter introduced in Equation 1, because flt(/) does not 
have exactly the same form as given in Equation 1. 

* Though their method also estimates the spectral index of 
the GW characteristic strain spectrum, we assume marginalization 
over this parameter in our discussion here. 
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the amplit ude. They also us e a method similar in 
concept to lJenet et al.l ()2005[ ) to attempt to detect 
the GW signal from binary SMBHs. 

First, the non-Gaussianity of the GW-induced ToA 
variations means that t he estimate of the intrinsic GW- 
induced variance of the Je net et al.l (|2005) statistic will 
be incorrect. This will affect estimates of the detection 
significance, particularly in the "strong signal regime", 
where the effects of GWs in the ToAs are large com- 
pared to all other noise proces ses. Second , the limit 
on the GW amplitude placed by Jen et et al.l ()2006D will 
be biased. The Jenet ct al. (2006) limit was placed by 
finding the GW amplitude for which 95% of simulated 
statistic values were above the measured value. The dis- 
tribution of their statistic derived using our simulations 
wo uld be different. R uling out a GW amplitude using 
the lJenet et al.l ()2006[ ) technique does not necessarily rule 
out a GW signal corresponding to our simulations with 
the same confidence. Finall y, our results indica t e that 
the likelihoods evaluated b y Ivan Haasteren et all (|2009f ) 
andi Pemorest et"aLl (|2012f) will also be biased, leading to 
a similar effect on GW amplitude constraints made using 
their methods. A definitive statement on the magnitude 
of the consequences of our simulations for current con- 
straints on the GW amplitude cannot be made, however, 
without fully considering the various PTA data analysis 
methods. 

6.2. Single GW source detection prospects 

We have established that in every frequency bin of the 
5-year datasets we consider, a few strong GW sources 
dominate the PSD of (5f . This means that we cannot 
consider the expected GW signal from binary SMBHs to 
form a background0 We briefiy consider the prospects 
for there being single sources of GWs that are de- 
tectable by PTAs. Various methods of detecting and 
characterising individual continuous sources of GWs with 
PTAs have re cently been presented ("Yardlev et al.ll2010l: 
Bovle fc PeJ I2OI0I: [Corb in & Cornish 2010: jLee erahl 
20111: iBabak fc Sesanal " l2012l: lEllis et al.ll20ia ). There 
are, however, few pre dictions for the exp ected numbers 
of detectable sources. iSesana et al.l (|2009f ) analysed sim- 
ilar binary SMBH population models to those consid- 
ered here to suggest that a 5-year ToA dataset would in- 
clude 5—10 single GW sources above the mean "stochas- 
tic background" level, mainly at GW frequencies greater 
than 10~^Hz. Their definition of a resolvable source as 
one which has a (mean) strain amplitude that is greater 
than the mean background level is conservative. This is 
because a PTA is capable of spatial, as well as frequency 
resolution. The background contribution per spatial res- 
olution element of a PTA will be less than the all-sky 
background level, resulting in a higher source amplitude 
to background ratio for a bright source located in the 
resolution element. 

The exact number of resolvable GW sources given a 
GW background level for PTAs depends on the partic- 

^ This result is analogou s to the c ase of the extragalactic back- 
ground light (Domi'ngucz ct al.|[20111 '). where the summed electro- 
magnetic radiation from AGN and star-forming galaxies is domi- 
nated by strong individual sources, behind which myriad further 
objects combine to form an apparently isotropic background too 
uniform to be resolved by current telescopes. 
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Figure 9. The average strain amplitudes of the three highest- 
amplitude binary SMBHs in frequency bins with < fc < 10 for 
each realization of the population. The strain amplitudes are ex- 
pressed as fractions of the mean summed amplitude of the remain- 
ing sources. We also show the 5th and 95th percentiles of the strain 
amplitudes, with their deviations from the means scaled down by 
a factor of 10. We made 300 realisations of the source population 
to produce this figure. As indicated in the Figure, squares (the 
solid line) depict the mean amplitudes of the strongest sources, 
circles (the dashed line) depict the mean amplitudes of the second 
strongest sources, and triangles (the dotted line) depict the mean 
amplitudes of the third strongest sources. 

ular search method. For example, IBovle fc PenI ()2010f ) 
suggest that a PTA composed of N pulsars could resolve 
up to 2N /7 sources per frequency bin. In Figure 9, we 
present a simple indication of the expected amplitudes of 
strong individual sources in the < A; < 11 spectral bins 
of our fiducial 5-year dataset. Using 300 Case R12 real- 
izations of the GW source population, we found the mean 
strain amplitudes in each spectral bin of the strongest 
three GW sources. We express these amplitudes as mul- 
tiples of the mean summed strain amplitude, /irost, of the 
remaining sources. The errors in the /irost values were not 
included in the error bars as they were very small. 

If we consider the sources besides the strongest three 
in a spectral bin to form a "background" 13 it is clear 
that, for spectral bins with fc > 2, three sources, on aver- 
age, produce the same total strain amplitude as the re- 
maining sources. Even for the fc = spectral bin, three 
sources are expected to produce more than half the total 
strain amplitude of the remaining sources. Indeed, the 
strongest source in the fc = spectral bin has an average 
strain amplitude that is ^ 0.35/ircst, which implies that 
a PTA which can resolve out two-thirds of the sky will 
detect equal contributions from the source and from the 
background. 

Blind searches for single GW sources with PTAs are 
therefore important. PTA data analysis methods that 
attempt to detect an isotropic component will not opti- 

This is by no means a rigorous definition of a background 
relative to the number of sources. The exact definition is dependent 
on the single source search method and the characteristics of the 
PTA. 
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mally recover the entirety of the GW signal from binary 
SMBHs, and could perhaps miss a large component of 
the signal for some realizations of the GW source pop- 
ulation. A careful consideration of the efficacy of GW 
background detection methods as compared to search 
methods for single sources, given the predicted source 
characteristics, is required. 

6.3. Limitations of our approach to modelling the GW 
signal from binary SMBHs 

A shortcoming of our approach towards modeling Sf, 
and indeed of all predictions for the GW signal from bi- 
nary SMBHs to date, is the assumption of circul ar or- 
bits for all binaries. Recent work (e.g., rSesanaj [20101 : 
IPreto et al.ll2011l: IKhan et al.ll201lD suggests that binary 
SMBHs emitting GWs in the PTA frequency regime 
will have highly eccentr ic orbits. The cand idate bi- 
nary SMBH OJ287 re.g.. lValtonen et al.llMl is in fact 
modeled with an orbital eccentricity of ~0.7. The GW 
waveform of an eccentric binary radiating in the PTA 
band spans many frequencies, and does not follow the 
frequency-time relation of Equation 15, or the depen- 
dence of the GW strain amplitude on frequency of Equa- 
tion 14. Therefore, if most binary SMBHs radiating in 
the PTA band are eccentric, the predicted mean spectral 
slope of the characteristic strain spectrum (Equation 1) 
will change. We will also need to account for the bi- 
nary eccentricity distribution in our predictions of the 
statistics of GW-induced ToA variations. Other effects 
that require further investigation are the effects of gas 
and stars on binaries. We do not include these effects in 
our model because the understanding of their combined 
contributions towards specifying the binary SMBH pop- 
ulation is not sufficiently advanced. 

7. CONCLUSIONS 

We have used a sop histicated model for galaxy evolu- 
tion ijGuo et al.ll2011[ ) to predict the distribution of bi- 
nary SMBHs radiating GWs in the PTA frequency band. 
By drawing lists of GW sources from this distribution, 
we simulate the effects of GWs from binary SMBHs on 5- 
year pulsar ToA datasets. We compare these simulations 
(Case R12) with simulated pulsar datasets containing the 
effects of an equivalent-amplitude GW signal modeled as 
a random Gaussian process (Case H09). We estimate the 
PSDs of the simulated GW-induced ToA variation time- 
series, and the correlations between these time-series for 
different pulsars. We find that the distributions of the 
PSD estimators of the realizations of the GW-induced 
ToA variations are different between the cases in every 
frequency bin, although the mean estimated PSDs are 
the same in each case. While in Case R12 the estimated 
PSDs are concentrated at lower values than in Case H09, 
the Case R12 estimations extend to higher PSD values 
than the Case H09 estimations. We also find that the 
functional form of the Hcllings & Downs curve is recov- 
ered on average in both cases. The correlations between 
the GW-induced ToA variation time-series for different 
pulsars in Case R12 are, however, significantly more scat- 
tered about the expected values than in Case H09. We 
interpret our results in terms of the influence of strong 
individual GW sources on the ToAs in Case R12. 

We conclude the following: 



1. The effects of GWs from binary SMBHs on pulsar 
ToAs cannot be accurately modeled using existing 
methods, i.e, as a random Gaussian process. This 
is because a few GW sources dominate the PSD 
of the GW-induced ToA variations at all frequen- 
cies, with reducing numbers of sources contributing 
equivalent PSD fractions in higher frequency bins. 

2. Our results directly affect existing PTA data anal- 
ysis methods aimed at detecting or estimating the 
parameters of the GW signal from binary SMBHs. 
The projected detection significance will be biased. 

3. The prospects for single GW source detection are 
strong. Individual sources could potentially be re- 
solved in all GW-dominatcd frequency bins of a 
5-year dataset. 

We emphasize that future searches for GW signals from 
binary SMBHs in pulsar datasets need to be sensitive to 
both individual sources as well as a GW background. 
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